library(tidyverse)

f <- function(x, pos) filter(x, `Number of Protein Chains in Target (>1 implies a multichain complex)` == 1, !is.na(`Ligand SMILES`))

t25.chembl <-
  read_tsv_chunked('~/BindingDB_ChEMBL.tsv', callback = DataFrameCallback$new(f))

# for data with polymer:
# 'Number of Protein Chains in Target (>1 implies a multichain complex)' ==1
keycol <- c('Ki (nM)', 'Kd (nM)', 'IC50 (nM)', 'EC50 (nM)',
            'BindingDB Reactant_set_id', 'Ligand InChI', 'Ligand SMILES',
            'PubChem CID of Ligand',
            'UniProt (SwissProt) Primary ID of Target Chain',
            'BindingDB Target Chain Sequence')

t25.chembl <-
  read_tsv('~/BindingDB_ChEMBL.tsv', col_select = keycol)

t25.chembl |> colnames()

t25.chembl <- t25.chembl |>
  drop_na('Ligand InChI', 'Ligand SMILES',
          'PubChem CID of Ligand',
          'UniProt (SwissProt) Primary ID of Target Chain')

t25.chembl |>
  mutate(`IC50 (nM)` = str_remove(`IC50 (nM)`, '>') |> str_remove('<') |>
                  as.double()) |>
  filter(!is.na(`IC50 (nM)`), `IC50 (nM)` < 1e7) |>
  select(-c('Ki (nM)', 'Kd (nM)', 'EC50 (nM)')) |>
  write_tsv('~/bindingdb.chembl.tidy.ic50.tsv')

t25.patent <-
  read_tsv('~/BindingDB_Patents.tsv', col_select = all_of(keycol))

t25.patent <- t25.patent |>
  drop_na('Ligand InChI', 'Ligand SMILES',
          'PubChem CID of Ligand',
          'UniProt (SwissProt) Primary ID of Target Chain') |>
  mutate(`IC50 (nM)` = str_remove(`IC50 (nM)`, '>|<') |>
           as.double()) |>
  filter(!is.na(`IC50 (nM)`), `IC50 (nM)` < 1e7) |>
  select(-c('Ki (nM)', 'Kd (nM)', 'EC50 (nM)')) |>
  write_tsv('~/bindingdb.patent.tidy.ic50.tsv')

t25.pubchem <-
  read_tsv('~/BindingDB_PubChem.tsv', col_select = all_of(keycol))

t25.pubchem |>
  drop_na('Ligand InChI', 'Ligand SMILES',
          'PubChem CID of Ligand',
          'UniProt (SwissProt) Primary ID of Target Chain') |>
  mutate(across(all_of(c('IC50 (nM)', 'EC50 (nM)', 'Kd (nM)', 'Ki (nM)')),
                \(x)str_remove(x, '>|<') |> as.double())) |>
  summarise(valid.ic50 = sum(!is.na(`IC50 (nM)`) & `IC50 (nM)` < 1e7),
            valid.ec50 = sum(!is.na(`EC50 (nM)`) & `EC50 (nM)` < 1e7),
            valid.Ki = sum(!is.na(`Ki (nM)`) & `Ki (nM)` < 1e7),
            valid.Kd = sum(!is.na(`Kd (nM)`) & `Kd (nM)` < 1e7))

t25.pubchem <- t25.pubchem |>
  drop_na('Ligand InChI', 'Ligand SMILES',
          'PubChem CID of Ligand',
          'UniProt (SwissProt) Primary ID of Target Chain') |>
  mutate(`IC50 (nM)` = str_remove(`IC50 (nM)`, '>|<') |>
           as.double()) |>
  filter(!is.na(`IC50 (nM)`), `IC50 (nM)` < 1e7) |>
  select(-c('Ki (nM)', 'Kd (nM)', 'EC50 (nM)'))

read_tsv('~/BindingDB_PDSPKi.tsv', n_max = 10)

t25.pdspki <-
  read_tsv('~/BindingDB_PDSPKi.tsv', col_select = all_of(keycol))

t25.pdspki |>
  drop_na('Ligand InChI', 'Ligand SMILES',
          'PubChem CID of Ligand',
          'UniProt (SwissProt) Primary ID of Target Chain') |>
  mutate(across(all_of(c('IC50 (nM)', 'EC50 (nM)', 'Kd (nM)', 'Ki (nM)')),
                \(x)str_remove(x, '>|<') |> as.double())) |>
  summarise(valid.ic50 = sum(!is.na(`IC50 (nM)`) & `IC50 (nM)` < 1e7),
            valid.ec50 = sum(!is.na(`EC50 (nM)`) & `EC50 (nM)` < 1e7),
            valid.Ki = sum(!is.na(`Ki (nM)`) & `Ki (nM)` < 1e7),
            valid.Kd = sum(!is.na(`Kd (nM)`) & `Kd (nM)` < 1e7))

read_tsv('~/bindingdb.chembl.tidy.ic50.tsv') |>
  bind_rows(t25.patent, t25.pubchem) |>
  write_tsv('~/bindingdb.3src.tidy.ic50.tsv')

t25.pubchem |>
  mutate(naa = nchar(`BindingDB Target Chain Sequence`)) |>
  ggplot(aes(naa)) +
  geom_density()

# performance check from Kexin H -----
dpperf <- read_tsv('../projects/DeepPurpose/performance.tsv')

dpperf <- dpperf |> mutate(
  MSE = str_remove_all(MSE, '\\s'),
  MSE.err = str_extract(MSE, '0.\\d+(?=\\))') |> as.double(),
  MSE = str_remove(MSE, '\\(.+') |> as.double(),
  cnci = str_remove_all(`Concordance Index`, '\\s'),
  cnci.err = str_extract(cnci, '0.\\d+(?=\\))') |> as.double(),
  cnci = str_remove(cnci, '\\(.+') |> as.double())

g1 <- dpperf |>
  mutate(model = fct_reorder(model, MSE, .desc = T)) |>
  ggplot(aes(x = MSE, y = model)) +
  geom_col(fill = 'salmon') +
  geom_errorbarh(aes(xmin = MSE-MSE.err, xmax = MSE+MSE.err), height = .3) +
  facet_grid(~dataset, scales = 'free_x') +
  theme_bw()

g2 <- dpperf |>
  mutate(model = fct_reorder(model, cnci)) |>
  ggplot(aes(x = cnci, y = model, fill = dataset)) +
  geom_col(fill = 'cyan3') +
  geom_errorbarh(aes(xmin = cnci-cnci.err, xmax = cnci+cnci.err), height = .3) +
  facet_grid(~dataset, scales = 'free_x') +
  labs(x = 'Concordance Index') +
  theme_bw()

library(patchwork)
g1 / g2

# DTC data -------
dtc <- read_delim('~/DTC_data.csv.gz')
dtc.prob <- problems(dtc)
dtc[183567, ]

dtc.prob |> count(col, expected, actual)

dtc |> tail()


panp <- read_delim('virtual_screen_v3/hpcc.adgpu.log', delim = ' ',
                   col_names = c('time', 'file'))
panp <- panp |>
  mutate(time = str_replace(time, '\\+', ' ') |> as_datetime(),
         id = str_extract(file, '\\d+') |> as.integer()) |>
  arrange(time) 

panp |>
  ggplot(aes(time, id)) +
  geom_point() +
  geom_path() +
  scale_x_datetime(date_breaks = '12 hours')

panp$id |> max()

tibble(i = 1:127) |>
  filter(!(i %in% panp$id))
